###
library(tidyverse)
library(modelsummary)
library(gt)
library(brms)
dat <- readRDS("../working/selection_data.rds")

###
dat$outcome <- as.numeric(dat$Funded == "Yes")
dat <- subset(dat, !is.na(ConMaj.allm))
mod_reward <- glm(outcome ~ log(Pop) + Rank + ConWinner1.allm + Region,
                  family = binomial,
                  data = dat)

mod_marg <- glm(outcome ~ log(Pop) + Rank + abs(ConMaj.allm) + Region,
               family = binomial,
               data = dat)

mod_comb <- glm(outcome ~ log(Pop) + Rank + ConWinner1.allm + abs(ConMaj.allm) + Region,
                family = binomial,
                data = dat)

AIC(mod_reward)
AIC(mod_marg)
AIC(mod_comb)

### Calculate Akaike weights
aic_vec <- c(AIC(mod_reward),
             AIC(mod_marg),
             AIC(mod_comb))

ev <- exp(-0.5 * aic_vec)
wts <- ev / sum(ev)

BIC(mod_reward)
BIC(mod_marg)
BIC(mod_comb)

BF_BIC = exp((BIC(mod_reward) - BIC(mod_marg))/2)  # From BICs to Bayes factor
BF_BIC  #

###
cm <- c("log(Pop)" = "log(Population)",
        "Rank" = "Rank within region",
        "ConWinner1.allmTRUE" = "In Conservative-held seat",
        "abs(ConMaj.allm)" = "Gap between Conservatives and nearest party")

gt_tab <- msummary(list("Core" = mod_reward, "Marginal" = mod_marg, "Combined" = mod_comb),
         coef_map = cm,
         output = "gt")

text <- gt_tab %>%
    tab_header(
        title = "Logistic regression model of town selection",
        subtitle = "Regional fixed effects omitted"
    ) %>% 
    as_latex()

text <- sub("ption*", "ption", text, fixed = TRUE)
writeLines(text, con = "../article/inserts/regs.tex")


mod_core <- brm(outcome ~ log(Pop) + Rank + ConWinner1.allm + Region,
                family = bernoulli,
                cores = 4,
                save_all_pars = TRUE,
                  data = dat)

mod_marg <- brm(outcome ~ log(Pop) + Rank + abs(ConMaj.allm) + Region,
                family = bernoulli,
                cores = 4,
                save_all_pars = TRUE,
                data = dat)


mod_comb <- brm(outcome ~ log(Pop) + Rank + abs(ConMaj.allm) + ConWinner1.allm + Region,
                family = bernoulli,
                cores = 4,
                save_all_pars = TRUE,
                data = dat)

bayes_factor(mod_marg, mod_core)
bayes_factor(mod_comb, mod_marg)
